Simulating the dynamics of relativistic stars via a light-cone approach 
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We present new numerical algorithms for the coupled Einstein-perfect fluid system in axisymmetry. 
Our framework uses a foliation based on a family of light cones, emanating from a regular center, 
and terminating at future null infinity. This coordinate system is well adapted to the study of the 
dynamical spacetimes associated with isolated relativistic compact objects such as neutron stars. 
In particular, the approach allows the unambiguous extraction of gravitational waves at future 
null infinity and avoids spurious outer boundary reflections. The code can accurately maintain 
long-term stability of polytropic equilibrium models of relativistic stars. We demonstrate global 
energy conservation in a strongly perturbed neutron star spacetime, for which the total energy 
radiated away by gravitational waves corresponds to a significant fraction of the Bondi mass. As 
a first application we present results in the study of pulsations of axisymmetric relativistic stars, 
extracting the frequencies of the different fluid modes in fully relativistic evolutions of the Einstein- 
perfect fluid system and making a first comparison between the gravitational news function and the 
predicted wave using the approximations of the quadrupole formula. 

PACS numbers: 04.25.Dm, 04.40.-b, 95.30.Lz, 04.40.Dg 



I. INTRODUCTION 

The dynamics of relativistic compact objects is of in- 
terest both to the astrophysical and relativity communi- 
ties. These objects represent on the one hand systems 
that may be active in stellar collapse, central engines for 
Gamma Ray Bursts, and are potentially strong sources of 
gravitational radiation. On the other hand, they repre- 
sent non-trivial solutions to the Einstein equations, pos- 
sibly exhibiting directly the predictions of that theory for 
cosmic censorship and black hole formation. 

Despite the obvious interest in building a reliable 
spacetime description of such systems, the process has 
been marred with difficulties |^]. One of the central 
problems concerns the formulation and solution of the 
initial value problem for the Einstein equations, a re- 
search effort that is still largely ongoing. The numeri- 
cal relativist must overcome a daunting list of problems, 
including a long term stable formulation of the initial 
value problem, the preservation of the constraints, the 
consistent imposition of boundary conditions and the re- 
liable extraction of physical information (including grav- 
itational waveforms). 

We propose that for a wide and interesting class of 
spacetimes most of the above concerns can be success- 
fully dealt with by the use of a characteristic initial value 
formulation for the Einstein equations. In particular we 
argue that the study of the non-linear dynamics of iso- 
lated relativistic stars is optimally performed within this 
framework and we proceed to illustrate the technical pro- 
cedures that support this claim. 

The incorporation of matter fields in the characteris- 
tic formulation of the Einstein equation was considered 
as early as 1983 || but the successful integration of the 



coupled system had to wait for the development of sta- 
ble algorithms for the vacuum Einstein problem. One 
dimensional schemes were developed in [Q. Algorithms 
for axisymmetric spacetimes, notably including a regu- 
lar origin, were presented in Recently using similar 
techniques, axisymmetric vacuum black hole spacetimes 
have been evolved §J. Techniques for extending finite 
difference algorithms to 3D were presented in [frl . Three 
dimensional codes excluding the origin of the light cones 
(hence presently unsuitable for stellar collapse studies) 
were presented in J§| . For an alternative approach see || . 

With reliable algorithms for the vacuum Einstein equa- 
tions available, a new line of research for the incorpora- 
tion of relativistic hydrod yna mics into numerical relativ- 
ity was initiated recently (lOj . This approach brings into 
the considerations the modern machinery from Computa- 
tional Fluid Mechanics. In this procedure, the evolution 
equations for the matter fields are solved using relativistic 
high-resolution shock-capturing (HRSC) schemes Jio]. [ll]] 
based upon (exact or approximate) Riemann solvers. A 
general formalism has been developed, which is now be- 
ing systematically applied to problems of increasing com- 
plexity. Applications in spherical symmetry have already 
been presented in the literature: investigations of accret- 
ing dynamic black holes can be found in |h], |l2) . Studies 
of the gravitational collapse of supermassive stars are dis- 
cussed in [^3| and studies of the interaction of scalar fields 
with relativistic stars are presented in jl4| . We note that 
there has been already a proof-of-principle demonstration 
of the inclusion of matter fields in three dimensions jlB) . 

In this paper we apply the general framework laid 
out in p"c| ] in the context of axisymmetric neutron star 
spacetimes. The code we have developed has been built 
upon previous work by Gomez, Papadopoulos and Wini- 
cour B, who constructed an axisymmetric character- 
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istic vacuum code. The numerical implementation of 
the field equations of general relativity is based on the 
light cone formalism of Bondi Jl6| , |Tt|] and Tamburino- 
Winicour flq] . 

The broad target of this project is the investigation 
of relativistic stellar dynamics (including collapse) using 
numerical relativity. A prerequisite for such studies is 
the development of very accurate and long-term stable 
general relativistic codes. It would seem that the fea- 
ture list of the characteristic approach makes it ideal for 
such studies. One serious bottleneck of the approach is 
the breakdown of a lightlikc coordinate system in the 
emergence of light caustics. Interestingly though, due 
to its quasi-spherical nature, no matter how agitated a 
relativistic core, it is unlikely to focus the light cones 
emanating from its interior. 

The light cone approach has a number of advantages 
compared to spacelike foliations: i) It is physically mo- 
tivated; the light cones offer a simple and unambiguous 
physical gauge on which to base the numerical spacetime 
grid, ii) It is unconstrained; the evolved variables capture 
rather directly the true degrees of freedom of the gravita- 
tional field, iii) It is very efficient; even in 3D, there are 
but two partial differential equations to solve, along with 
a set of radial integrations along the light cones, iv) It 
allows for well defined compactification of the do- 
main, which leads to perfect outer boundary conditions, 
v) Finally, and perhaps most importantly, the above the- 
oretical advantages have been shown in a series of works 
to translate to remarkably numerically robust and stable 
codes (see e.g. |2(|). For a recent review of the approach 
the reader is referred to 21 1. 

An important focus of research in numerical relativ- 
ity has been the investigation of dynamical spacctimes 
of relativistic stars. While the pioneer investigations 
were mostly concerned with the study of spherically- 
symmetric gravitational collapse scenarios the effort is 
nowadays expanding to the study of the merging of com- 
pact neutron star binaries. Advances in the numerical 
formulation of the equations as well as in computational 
power have made possible to transform the numerical 
evolution of neutron stars in general relativity into an 
active field of research || ||, |J, ||, |ej. 

The paper is organized as follows: In Sec. II we de- 
scribe our mathematical framework for characteristic nu- 
merical relativity in axisymmetry. The next section de- 
scribes our numerical implementation in some detail. In 
Sec. |^we present tests to assess the validity of the differ- 
ent regimes of our numerical implementation. In Sec. ^ 
we apply the current code in studies of stellar pulsations. 
The last section concludes with a short summary. 



V a T Qb = 
V a {pu a ) = 0, 



(2) 
(3) 



with the latter two equations being the local conservation 
laws of stress-energy and current density. The energy- 
momentum tensor T ab has the form 



T ab = phu a u b + pg a b- 



(4) 



In this expression p denotes the rest mass density, h = 
1 + e + 2 is the specific enthalpy, e is the specific internal 
energy and p is the pressure of the fluid. The four-vector 
u a , the 4- velocity of the fluid, fulfills the normalization 
condition g a bU a u b — — 1. Using geometrized units (c = 
G = 1) the coupling constant is k = 8ir. In order to close 
the system of fluid equations, an equation of state (EoS) 
has also to be prescribed, p — p(p, e). 



A. The Einstein equations for the Bondi metric 



As a geometric and coordinate framework we will use 
the Bondi (radiative) metric Jll| to describe the space- 
time 



ds 2 



V 



e 2p _ u2 r 2 e 2n du 2 _ 2e ^dudr 



2Ur 2 e 2 ~<dud9 + r 2 {e 2l dQ 2 + e' 2 ' 1 sin 2 6>d</> 2 ), (5) 



with coordinates (x 3 ) = (u,r,9,(f>), where u is 

a null coordinate labeling outgoing light-cones, r is the 
radial coordinate whose level surfaces (two-spheres) have 
area 4irr 2 , and and <fi are angular coordinates with <fi 
being a Killing coordinate. The metric functions V, U, (3 
and 7 depend on the coordinates u, r and 9. We choose 
the origin of the coordinate system r = to lie on the 
axis of our axisymmetric stellar configurations. With the 
above assumptions, the gravitational field equations read 
explicitly 



II. MATHEMATICAL FRAMEWORK R ab = K ( ph ( u a u b + \g ab ) - pg ab ) , (6) 
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We will work with the coupled system of Einstein and 
relativistic perfect fluid equations 



where the relevant components of the Ricci tensor R a b 
G ab = KT ab , (1) are 
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~ A Rrr = Ar - ^{l ,rf , (7) 

2r 2 R r0 = (rV^-^Xr - 2r 2 (/3 r9 - 7 , r9 + 2 7 . r7 , e - -/3 e - 2 7 , r cot0) , (8) 

r 

-rVV*i^ ^ 2^ + irV(-«(^-^, p9 -4r^-^, P cot fl -4r^ cot. (9) 

+2e 2(/3-7){_! _ (3 7 _ e _ cot6 _ 7 + ^ + (/3 e )2 + 27 e(7 e _ p ^ i 

- r*(Mg++Rto = 2r(r 7 ),„ r + (1 - r 7r )V> - (r 7>rr + 7 , r )V - r(l - r 7 , r )£/ e - r 2 (cot - %e )U r 
+e 2(/3-7)(_! _ (3 7 e _ 2/ 3 ifl ) cot - 7, ee + 2 7>e ( 7 , e - t0 )) 

+r U(2Ty >r e + 2 7 ,e + r 7 , r cot - 3 cot 0) . (10) 



In Eq. @ A, _B denote the angular coordinates, A, B = 
2,3. As usual a comma is used to denote a partial deriva- 
tive. 

The Einstein equations decompose into hypersurface 
equations, evolution equations and conservation laws. 
The hypersurface equations, Eqs. (@)-(H), form a hier- 
archical set for P <r , U. r and V. r . The evolution equation 
is an expression for (r 7 ), ur given by Eq. (|Io|). The light- 
cone problem is formulated in the region of spacetime 
between a timelike worldtube T, which in our case is lo- 
cated at the origin of the radial coordinate r = 0, and 
future null infinity J + . Initial data 7 is prescribed on an 
initial light cone u = in this domain. Boundary data 
for P, U, V and 7 is also required on T. 

As shown in the original paper of Bondi jl6| , the con- 
tracted Bianchi identities for the vacuum field equations 
enforce all other Ricci tensor components to vanish, if 
they vanish on a worldline. In the same way one can 
show, that the contracted Bianchi identities for the mat- 
ter system enforce all other components of the Einstein 
equation (||) (see (|]). 

The reformulation of the above form of the equations 
for numerical integrations follows the work of Ref. ||. 
First, in order to be able to compactify the entire space- 
time and to better resolve interesting parts of the nu- 
merical integration (e.g. a stellar configuration centered 
at r — 0) we allow, starting from a radial coordinate 
x G [0, 1], for a coordinate transformation of the radial 



coordinate x — > r(x), with an associated derivative 

dx/dr=f(x). (11) 

This transformation generalizes the results of @], where 
the fixed grid r(x) — x/(l — x) was used. Furthermore, 
in order to eliminate singular terms at the poles we use 
the new coordinate |27l 



y = — cos a 

and we introduce the new variables 

V-r 



u = 



U 
sin0' 

7 

sin 2 6 ' 



(12) 

(13) 
(14) 

(15) 



The metric, hence, takes the form 



ds 2 



V 



-e 2/3 + C/Ve 27 du 2 - 2f- 2 e w dudx 



2Ur 2 e 2l dudy + r 2 (e 27 sin~ 2 8dy 2 

2 ). (16) 



-e~ 27 sin 2 



The hypersurface equations (@)-(j|) thus read 



(W 



2 p 2( 7 y-/3) 



u x 



2/V(7,x) 2 + i / 2 i? M , 



2r \ /3 tXy - -^pP,y + A Vl,x + J/[27,s(j/7,« - 2y 7 ) - j <xy ] \ + 2r R xy , 



(17) 
(18) 



r 2 f 2 S, x + 2rS = -1 - AryU - r 2 fyU, x + 2ryU y + y (jf 2 U, xy ~ ^/V^-«(*7 X ) 2 ) 



- 12 7 - 2y/3, y + y[10 7 + 8yj, y + 8 7 2 + 4y 7/ 3, + p. yy + (A,) 2 ] (19) 



f[%f + 2%P, y + \ yy + 8yn y ] + 2y 3 7f„} - ^g AB R AB 
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We have used here the notation y = 1 — y 2 and and <? a v a V b ?/; = —e~ 2l3 H, (22) 

O.j, to denote the partial derivatives with respect to x 
and y, in contrast to the partial derivatives with respect 

to r and 6. Note, that the Ricci tensor components are where the ~ quantities are induced by the 2-metric 
those in x, y-coordinates. 

Due to this choice of variables, the Einstein equations y 

are non-singular on the polar axis, where y = ±1. Note, da 2 = e 2/3 du 2 — 2e 2/3 du dr, (23) 

that the y-component of the four-velocity fulfills 



r 



(20) i.e. explicitly 
smtr 



which is in analogy to Eq. (|14|). 



V , 



The evolution equation for 7 is written in the form of g a V a Vfc?/; = — e 1 {2ip ur — ( — Vv),r)i (24) 



a wave equation for the quantity tp 

ip = ry, (21) and 



r 



~(S + rf 2 S, x ) 7 + l -e 2 ^-^ [fi m + (f3 >y ) 2 ) + 4jUy 
+2rf 2 ^U, x y + 6rf 2 %Uy - y (rfl, y ll, x + 2rf\ xy U + 2y lV Uj 

+ ]_ e 2(()-iy) KphUyUy ( 25) 



B. The relativistic perfect fluid equations 

Whereas Eq. (§) is a strict conservation law, Eq. (12j) 
involves source terms when writing the covariant deriva- 
tives in terms of partial derivatives. In the presence of a 
Killing vector field it can be recast as a conservation law. 
Following p8[ , the number of source terms in Eq. (|J) is 
minimized using the equivalent form 

V a T a b = 0. (26) 

In our hydrodynamics code we use the form given by 
Eq. (^) in order to set up the evolution equation for the 
radial momentum. This is motivated by stability consid- 
erations when evolving spherical neutron star models (see 
below). However, to set up the evolution equation for the 
polar component of the momentum we use Eq. ( p6| ) . This 
form of the conservation law eliminates singular behavior 
of the y-component of the velocity at the polar axis. 

After introducing the definitions U° = ^—gT 00 , 
U x = V=gT 0x , U y = ^/=gT° y and U 4 = ^=gpu°, 
the fluid equations can be cast into a first-order flux- 
conservative, hyperbolic system for the state- vector U = 

(U ,U* : Uy,U 4 ) 

doU + djF j0 = S° , (27) 
d U x + djF jx = S x , (28) 



d U y + djF> y = S y , (29) 

d U 4 + djF^ = S 4 . (30) 
The flux vectors are defined as 

F j0 = J~g T j0 , (31) 

F rx = y/=g T 1X , (32) 

F\ = V=9~T{, (33) 

F ji = V^gpu* , (34) 
and the corresponding source terms read 

S a = g ab S b + J=g~T c b d c g ab , (35) 

&a = V~ 9 T J-ab 

= -^-^-phu c u d (g cd ), a +p(V Z lj),a , (36) 

S 4 = , (37) 



In the above expressions \J—g is square root of the four 
dimensional metric determinant and T bc are the Christof- 
fel symbols. 

C. Gravitational waves at null infinity 

Using a compactified coordinate x E [0,1], 
lim x _>i r(x) — 00, we have future null infinity J~ + on our 
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grid, where we can unambiguously extract waveforms. 
Hence our approach does not suffer from the problem 
most numerical relativity codes have to deal with, i.e. 
extracting approximate gravitational waveforms at a fi- 
nite distance (for a comparison, see 29, 30]). Performing 



a power series expansion of 7 around future null infinity 
in inverse powers of the radial coordinate r 



7 = K + - + 0(r~ 2 ), 



(38) 



the hypersurface equations yield an expansion 



/3 = H + 0(r 



(39) 



U = L + Oir- 1 ), (40) 
V = r 2 (Lsm6), e /sm6 + re 2(H - K) [l +K >ee 

+2{H fi sin ff) >6 / sin 9 + 3K,g cot 9 + A{H fi f 
-\R fi K fi - 2(K,e) 2 ] - 2e 2H M + O^ 1 ), (41) 



where M = M(9) denotes the Bondi mass aspect. As 
its straightforward extraction at J + will be spoilt by the 
leading, diverging terms in Eq. (^l]), we follow the pro- 
cedure proposed in |3l) to determine the Bondi mass of 
a numerical spacetime. We numerically solve the hyper- 
surface equations for the new metric variables 



2r = (l-y 2 r 1 / 2 r 3 e 2 ^ m U, r + 2rf3,y 
2/i = -V + r 2 [(l-y 2 ) 1 / 2 U}, y +r 3 e 2 ^ 



-^(I-jT^KI 



— (l-y 2 )e- 2 ^ 



+ e 2 ^[(l-y 2 )e 



■ y 



(42) 
(43) 



r 



(see ||l| for more details) . The Bondi mass M B can then 
be readily computed as 

M B = — [ uj^e- 211 uL =1 sin<9 dO dd>, (44) 
An J 

where u> denotes the conformal factor relating the two- 
geometry 



ds 2 = e 2K d9 2 + sin 2 9e~ 2K dd> 2 



to the two-geometry of a unit sphere 



ds 2 B = d9 2 B + sin 2 9 B d(f) 2 B . 



(45) 



(46) 



The total energy emitted by gravitational waves during 
the time interval [iijU + du] in angular directions [9, 9 

dO] x [(f>,(f> + d(f>] is 



dE = —t—uj ~e 

l07T 



1 211 J , (sin6> c 2 L), g 



+e- 2K Lusin9[ 



(e 2 "w), 



sin v c 
2 

sin0d6>d0du.(47) 



The quantities K, c, -ff and L are read off from the metric 
variables at J + , e.g. c = — {r 2 f 2 ^)\ x= i (our coordinate 
transformations r = r(x) fulfill the requirement that r 2 / 2 
is finite). 

For the extraction of waveforms seen from a distant 
inertial observer, we have to transform our coordinate 
system to a Bondi coordinate system. Following || the 
Bondi coordinate time u B is related to the retarded time 
u as 



1 



du = — e du B , 



(48) 



whereas the angular Bondi coordinate y B = —cos(9 B ) 
can be calculated from 



dy = ~^dy 



With the definition of the news function 

1 e- 2H [ (sin 9 c 2 L), e 

2 uj z 



-IK ■ a A e ^)fi 



(49) 



(50) 



and after integrating over the angles <fi = <fi B , one recovers 
the expression for the total energy radiated according to 
Bondi 111 



dE = ^N 2 dy B du B . 



(51) 



This relation allows us to check global energy conserva- 
tion 



ec := M B (u) - M B (u = 0) + 



[ [ dE = 0. (52) 
Jo J -1 



For the calculation of the Bondi mass, Eq. (|4j), as 
well as for the calculation of the Bondi news, Eq. j5C|), 
we have to determine the conformal factor oj 



ds 2 B = ui 2 ds 2 . 



(53) 



It can be shown that the coordinate y B = — cos(9 B ) reads 

J2K 



y B (y) = tanh - 



1 f v e 2K -I ,„ 1 f v e 2K - 1 



-^-dy + - 



1-y* " 2J 1 1-y 



^rdy 



i-y 



-^dy 



(54) 
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The choice of the integration constants ensures regularity 
of db, i-e. Um y ->±iyB = ±1 and for spacetimes with 
equatorial plane symmetry, ys is symmetric as well. The 
conformal factor can be written as 



2c 



K 



(1 + y)e A + (1 - y)e- A ' 



(55) 



where 



A(y) 



V C 2K 



-=v~dy 



2 7_i i-y 2 2 A 1-y 



y „zk 



^rdy. (56) 



The regularity of the conformal factor can be directly 
seen if we write Eq. ( [slf ) as 

rV /•! fV /■! 

A(y) = dy dae 2aK K + dy dae 2aK K, 



J-i Jo Ji Jo 

where using Eq. ( |l5[ ) we have defined 

K 



K 



sm 



(57) 
(58) 



which makes explicit use of the full spectral decomposi- 
tion of the system. For our particular formulation of the 
hydrodynamic equations such characteristic information 
was presented in [^o| . We note that unlike in previous 
work JT^I , we have now included the metric determinant 
in the definition of the conserved quantities U. Never- 
theless, as was explicitly demonstrated in |2S[] , the eigen- 
values needed for the computation of the local Riemann 
problems remain unchanged. 

In more precise terms the hydrodynamics solver of our 
code uses a second order Godunov-type algorithm, based 
on piecewise linear reconstruction procedures at each cell- 
interface J32] , |33| and for simplicity the HLLE approxi- 
mate Riemann solver |34|, ^5) . For the time update we use 
the second order Runge-Kutta algorithm derived in |36| . 
General information on HRSC schemes in relativistic hy- 
drodynamics can be found, e.g. in [ pT[ |3?| and references 
therein. 

We use the procedure described in Q to explicitly 
recover the physical (primitive) variables for the perfect 
fluid EoS p = (T — l)pe, T being the adiabatic index of the 
fluid. We solve for the primitive variables e, p, u x and u y 
and then use the normalization condition to determine 



III. THE NUMERICAL IMPLEMENTATION 



We use an equidistant grid covering our numerical do- 
main (x,y) £ [0,1] x [—1,1] with grid spacings Ax = 
1/N X , Ay — l/Ny, where N x + 1 is the number of grid 
points in the radial direction and 2N y + 1 is the number of 
grid points in the angular direction (N y is the number of 
angular grid zones per hemisphere) . All variables are de- 
fined on the grid (u™ ,Xi,yj) = (nAu,iAx,jAy), except 
for the quantities U and r which are defined on a stag- 
gered grid (u n , x l+1/ 2,y j +i/ 2 ) = (nAu, (i + 1/2) Ax, (j + 
1/2) Ay). As shown by || the use of a staggered grid is 
necessary for stability. 



B. The metric evolution 

For the time update of the metric field 7, we solve 
the wave equation (22) with a so-called parallelogram al- 
gorithm, based on a parallelogram consisting of ingoing 
and outgoing characteristics jj, |5[ [38] . We closely follow 
the implementation of jj| to which the reader is referred 
for more details. In contrast to this work, however, we 
only use a single form of the wave equation in our nu- 
merical implementation, which is regular at J + . This 
regularization is accomplished by rewriting the parallel- 
ogram identity in terms of the quantity 



The fluid evolution 



(60) 



In our implementation of the fluid equations we closely 
follow the work of Papadopoulos and Font [jl(|. The 
hyperbolic mathematical character of the hydrodynamic 
equations allows for a solution procedure based on the 
computation of (local) Riemann problems at each cell- 
interface of the numerical grid. At cell i,j the state- 
vector U is updated in time (from u n to u n+1 ) using a 
conservative algorithm 



U" +1 



U" 



^(G ij+ i/2 - G ij _ 1/2 ) 
AuSjj , 



(59) 



where the numerical fluxes, F and G, are evaluated at the 
cell interfaces according to some particular flux-formula 



Explicitly, for a parallelogram S with left upper corner 
P, right upper corner Q, left lower corner R and right 
lower corner S (see Fig. m the algorithm reduces to 



- Au r Q f\ Q 


H c 








f\ P ( p 4 


Au 


rpf\p 


He) 




/Is 1 6 4 


Au 


rsf\s 


H c ) 




f\ R { R 4 


Au 


rnf\n 


He). 


(61) 



In the above expression H c denotes the source term of 
Eq. ( p5| ) at the center of the parallelogram approximated 
to second order accuracy. The regularity of this expres- 
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R 



A. 



i+1 



FIG. 1: Parallelogram PQRS consisting of ingoing and out- 
going light cones. We have displayed the outgoing light cones 
at two different coordinate values u, A and B, foliating our 
spacetime. For simplicity, all null lines have angles 45°. In ad- 
dition to the parallelogram we have marked radial grid points 
at i and i + 



sion can be easily seen for our compactified grids 

x 



as 



and 



r{x) 



rf 



1 



b > 1, 



>/! + (& -l)a 



(62) 



(63) 



f\Q 
/Is 



1-XQ S?=o4 Vl + ft-W 

E;=o4 ,/l + (&-!)*» 



(64) 



All the terms on the right hand side of Eq. (j64| ) are ex- 
plicitly regular for x — > 1 (r — * oo), except for the first 
factor which is regular according to the discussion in || . 

We discretize the hypersurface equation ([l7]) (and sim- 
ilarly Eq. ©) as 



(65) 



where Tip denotes the right hand side of Eq. (|l7|). To 
solve the hypersuface equation (18), we discretize the al- 
ternative equation 



2 f 2 



2xf [r*ru x ) x4 +r z f 



i 



3„2(/3-y7) 



H w , (66) 



where the right hand side of Eq. ( Eq ) has been denoted by 
2r 2 Tiu. The derivative d/dx 4 — ^^d/dx was introduced 
to ensure regularity at the origin. 

In the following, we describe the order of the time up- 
date from light-cone A at time u to light-cone B at time 
u + Au (see Fig. 0). Let us assume that we know the 
primitive and conserved fluid variables, and the metric 
quantities 7, (3, U and S on the light cone A. In a first 
step we globally determine the conserved fluid variables 
on B. For the metric update, in contrast, we march from 
the origin to the exterior of the light-cone B. Having 
previously obtained the variables up to grid point i (ei- 
ther from the specific boundary treatment at the origin 
or during the marching process), we first determine 7 at 
grid point i + 1. In a second step we solve for /3, U and 
S at i + 1 in that particular order, recovering the prim- 
itive variables with the metric thus obtained. As the 
hypersurface integration for the metric depends on the 
primitive variables at the grid point to be determined, 
we iterate the hypersurface and recovery algorithm until 
convergence. 

In order to obtain stability in our explicit algorithms 
when solving the fluid and metric equations we have to 
fulfill the Courant-Friedrichs-Levy condition - the numer- 
ical domain of dependence must include the analytical 
domain of dependence. This limits the maximal time 
step allowed in each time update. Calculating the char- 
acteristic speeds for the fluid system, the fluid update 
sets a limit on the time step as 



Ait < min(ciAa;, CaAy), 



(67) 



where c\ and C2 are constants and the minimum is calcu- 
lated for the entire fluid grid. For the metric update, it 
can be shown that the evolution near the origin sets the 
stricter theoretical limit 



Am = c 3 Ar(Ay) 2 



(68) 



with C3 = 0.5. In numerical experiments we found, how- 
ever, that with our coupled code C3 ~ 10, in good agree- 
ment with the result C3 = 8 reported injE7|. With this 
result the time step restriction from Eq. ( |68| ) is not much 
stronger than the time step restriction from Eq. ( |6^ ) , at 
least for the angular resolutions we can afford. In the 
simulations described in this paper we use time steps 
of 0.6 times the maximal time step consistent with the 
fluid evolution Eq. (|67]). With this procedure, it was not 
necessary to implement implicit methods for the metric 
update. 

As for the vacuum equations the origin of coordinates, 
where we assume that the coordinate system is a local 
Fermi system, needs special care. The main change in 
the falloff behavior of the metric variables due to the 
presence of the fluid is in (3. We impose a falloff behavior 
of the metric field 7 as 



7 = a r + 



b r d 



(69) 



where a and b are constants. At the origin the radial 
dependence of j3 is = 0(r 2 ), instead of j3 = (9(r 4 ) for 



8 



the vacuum case (from Eq. (|17]) ) . To ensure regularity of 
the fluid at the origin, we have to impose 0] 



u r = D + Ey + 0(r), 
u y = Er + 0(r 2 ), 



(70) 
(71) 



where D and E fulfill -D 2 + E 2 = -1. With the def- 
initions (the leading terms can be extracted from the 
hypersuface equations) 



where Y — Ve~ 2/3 . Moreover, for the study of oscilla- 
tions of these spherical stellar models discussed below, 
we give some additional fluid perturbation. In order to 
obtain consistent initial data, we have to solve the hy- 
persurface equations (ff|)-(||) imposing the normalization 
condition for the fluid. As we have presented results for 
a similar code elsewhere |14j , we keep the presentation of 
this section short. 



-nph{D + Ey) 2 + F{y)r 3 



3 1 

U = Ay(ar + -br 2 ) + -nph{D + Ey)Er 
5 2 

+C(y)r 2 , (73) 
KphUyUy = nphE 2 r 2 + G(y)r 3 , (74) 

the quadratic terms in r in the wave equation reduce to 
an equation for a 



\F.yy + ^G. 



(In the absence of matter this reduces to a tU = 
see ||). We extract a, b, C, F and G at the old time 
slice and then solve Eq. ( f75| ) to obtain a at the new time 
slice. Inserting this value into the leading order in r of 
Eq. (|6|) and Eq. ((73|), we calculate 7 and U at the two 
first grid points, which then allows us to start the march- 
ing algorithms described above for the metric update. 

In our numerical implementation we use a standard 
second order integration for the double integrals in 
Eq. ( p^ ) and we determine the conformal factor accord- 
ing to Eq. ( |55|) afterwards. Furthermore, to solve the 
hypersurface equations for the auxiliary variables r and 
p (n is needed for the Bondi mass, Eq. (^)), we also use 
a second order discretization. By simply reading off the 
metric values at J + , it is finally possible to determine 
the Bondi mass and news. 



(72) 1- Stable relativistic stars 

First we study equilibrium models of relativistic stars. 
In order to obtain initial data we solve the TOV equa- 
tions ([76|)-(|77]) for a polytropic EoS p — Kp r , K being 
the polytropic constant. With complete initial data at 
hand we evolve stable equilibrium configurations in time. 
Both, for finite grids which only cover the star, as well as 
for compactified grids, where we cover the star and its en- 
tire exterior spacetime, we are able to maintain the initial 

(75) 

equilibrium profiles of the star during a time-dependent 
simulation, for times much longer than the light-crossing 
time. Due to the discretization error, the stars are exited 
to oscillate in their radial modes of pulsations (see p3] 
for more details). Deviations from equilibrium converge 
to zero with increasing resolution with a second order 
convergence rate. 



IV. CODE TESTS 

In order to validate the accuracy of our numerical code 
we have performed various tests aimed to check the dif- 
ferent regimes of the implementation. 



A. Spherically symmetric tests 

In this section we describe tests of our code in dealing 
with spherically symmetric stellar models. In order to 
set up equilibrium models for relativistic stars we solve 
the Tolman-Oppcnheimer-Volkoff equations in outgoing 
null coordinates JlIJ 



P.r 



1 1 

27 ~ 2Y 
Y r = l + 87rr 2 (p- ph), 



(1 + 8irr 2 p) ph 



2. Migration of an unstable relativistic star 

Following |}9| we have checked the code on the dy- 
namical evolution of an unstable spherical star. In such 
a model the sign of the truncation error of the numerical 
scheme controls the fate of the evolution: the star may 
either expand or collapse. In our code this sign is such 
that the unstable star "migrates" to the stable branch of 
the sequence of equilibrium models. In such a situation, 
the rest-mass of the star has to be conserved throughout 
the evolution. Despite being an academic problem this 
simulation represents an important test of the accuracy 
and self-consistency of the code in a highly dynamical 
situation. 

As in (3^] we have constructed a n = 1 (r = 1 + 1/n = 
2), K = 100 polytropic star with mass M = 1.447 M Q 
and central rest-mass density p c — 8.0 x 10~ 3 (in units 
in which G = c = Mq = 1). As the radius of the star 
strongly increases during the evolution, we surround the 
star by a low density atmosphere. In order to avoid nu- 
merical problems in these zones we reset the fluid vari- 
ables to their original atmosphere values once they have 
fallen below a threshold value (see |40| for more details). 
This artificial resetting enables the star to expand and to 
contract. Despite the use of an atmosphere the energy 
conservation properties are well satisfied [H . 

(76) Fig. [2] shows the evolution of the central density up 
to a final time of ub = 10 ms. On a very short dy- 

(77) namical timescale the star rapidly expands and its cen- 
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0.01 



0.008 



0.004 



of spacetime. 




4 6 
Bondi time [ms] 

FIG. 2: Evolution of the central rest-mass density during 
the migration of an unstable relativistic star (n = 1,K — 
100, M = 1.447M 0) p c = 8.0 x 10~ 3 ;G = c = Af Q = 1) to 
a stable model with the same rest-mass. The central density 
of the (final) stable configuration is p c = 1.35 x 10 -3 . The 
evolution shows the expected behavior. Since we are using a 
polytropic EoS, the amplitude of the oscillations is essentially 
undamped for the evolution times shown. 



Exact vacuum solution SIMPLE 



Following the work of || the metric 



U 
V 



2 (1 + E) ' 

4S ' 
a 2 g cos 9 



(2aV ~aV + 1), 



is a solution of the vacuum field equations for 



g = rsm9, 



(78) 

(79) 

(80) 
(81) 



(82) 
(83) 



where a is a constant. Using this solution with a suitable 
L2-norm to measure deviations we have checked that our 
metric solver is second order convergent. 



tral rest-mass density drops well below its initial value, 
less than p c = 1.35 x 10 -3 , the central rest-mass density 
of the stable model of the same rest-mass. During the 
rapid decrease of the central density, the star acquires a 
large radial momentum. The star then enters a phase 
of large amplitude radial oscillations around the stable 
equilibrium model. As Fig. || shows the code is able to 
accurately recover (asymptotically) the expected values 
of the stable model. Furthermore, the displayed evolu- 
tion is completely similar to that obtained with an inde- 
pendent fully three-dimensional code in Cartesian coor- 
dinates |3j|. 

The evolution shown in Fig. ^ allows to study large am- 
plitude oscillations of relativistic stars, which could occur 
after a supernova core-collapse |4l| or after an accretion- 
induced collapse of a white dwarf. 



B. Tests beyond spherical symmetry 

To the best of our knowledge we do not know of any 
regular exact (analytic) solution of the Einstein equa- 
tions in axisymmetry with a non-vanishing perfect fluid 
matter field. However, there is an exact solution for the 
vacuum equations which was already used to check the 
vacuum code of |J]. As we have generalized the coordi- 
nate system, we use here the same solution to check our 
code. 

In addition, to test the overall implementation of the 
code we examine the global conservation properties which 
can be rigorously established due to the compactification 



2. Global energy conservation test 

In this section we focus on a global energy conservation 
test. Starting with the equilibrium model of a neutron 
star (see next section for more details), we use a strong 
gravitational wave to perturb the star, 



7 = 0.2 e 



-3(r-4)' 



(84) 



Such a large amplitude is not realistic, but we choose 
it to test our numerical implementation in the nonlinear 
regime. Fig. || shows a contour plot of the initial gravi- 
tational wave shear. In what follows we use a radial grid 
with r = Srr 9 —- 

Fig. |4| shows the deviation from global energy conser- 
vation (see Eq. (|52])) as a function of the grid resolution 
(circles) and the total energy radiated away in gravita- 
tional waves. The deviations from exact energy conser- 
vation converge to zero which represents a very severe 
global test for our numerical implementation. 

There are two additional consistency conditions, which 
relate the metric quantities on J + . These are |27| 



S-U.t 



7,u + ~e 27 sin ( 



,2~ 



U cot 
U 



sin 8 



= 0, 
= 0. 



(85) 
(86) 



Our code reproduces these conditions, the errors converg- 
ing to zero with a convergence rate of 1. Fig. || shows the 
deviation from zero for the first condition, thus checking 
the leading term in the falloff behavior of the quantity S 
at J+. 
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gravitational wave shear 



-0.20 -0,10 0.00 0.10 0.20 




FIG. 3: Contour plot of the initial gravitational wave shear 7, 
Eq. (|84|). The axes labels denote the value of the radial coor- 
dinate x along the equator (horizontal axis) and the symmetry 
axis (vertical axis), and they run from the origin of the coor- 
dinate system at x = to future null infinity J + at x = 1. 
The numerical domain comprises the half circle to the right 
of the vertical line at x = (the symmetry axis) , the left half 
is obtained from axisymmetry. The gravitational wave shear 
perturbation is located in a small ring at radius r « 4. 



The obtained first-order convergence rate can be ex- 
plained by the use of a total variation diminishing HRSC 
scheme for the fluid evolution, which, although it is 
second-order accurate in smooth, monotonous parts of 
the flow, reduces to first-order at local extrema, which 
are present in the interior of the numerical domain in 
this test (see [22 for alternative essentially nonos dilatory 
schemes). Tests including propagation and scattering off 
the origin of pure vacuum gravitational fields yielded the 
expected second-order convergence. 



V. PULSATIONS OF RELATIVISTIC STARS 

Oscillations of neutron stars, although believed to emit 
only weak gravitational waves (at least in the absence of 
instabilities) are interesting sources for future gravita- 
tional wave detectors (for a recent review see [Egj ). 

In this section we apply our code to the study of ax- 
isymmetric pulsations of relativistic stars. We focus our 
analysis on one specific stellar model. We choose a rel- 
ativistic polytrope p = Kp 1+ ~ with polytropic index 
n = 1, polytropic constants K = 100 and central density 



1.28 x 10 3 (in units in which c = M G 



Pc 

This model, which has a total mass of M 
already been used in previous work 22 , 23|, [39| , which al 



= G = 1). 
1.4M Q has 



4.0e-05 



total energy radiated 
; global energy conservation ec 




0.002 0.0025 
resolution 1/Nx 



FIG. 4: Global energy conservation test for a neutron star 
and a strong gravitational wave. Plotted are the deviation 
from global energy conservation ec (circles, see Eq. (|j^)) and 
the total energy emitted by gravitational waves as a function 
of the grid resolution. The final integration time is ub = 
2 10~ 8 s, the total number of angular grid points 2N y = 
0.1N X . 



0.0015 




angular coordinate 



FIG. 5: Consistency check for the global norm Eq. (|8fj) at J + . 
Plotted are deviations from zero as a function of the angular 
coordinate y and for the same resolutions used in Fig. ^| The 
errors decrease with resolution, the convergence rate is one. 



lows us to compare our results for the radial frequencies 
and fixed background evolutions in axisymmetry. 



A. The perturbations 

In order to excite the radial oscillation modes of the 
star (I = 0) we perturb the equilibrium configuration 
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anqular velocity 
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FIG. 6: Contour plot of the angular velocity perturbation u y , 
Eq. (|qo|) - The axes labels denote the value of the radial coordi- 
nate x along the equator (horizontal axis) and the symmetry 
axis (vertical axis), and run from the origin of the coordi- 
nate system at x = to future null infinity J + at x = 1. 
The star corresponds to the inner circle, where the velocity 
perturbations are non-zero, with the north pole above. 



TABLE I: Mode frequencies obtained in the Cowling approx- 
imation for a relativistic polytrope with K = 100, n — 1 
and central density p c = 1.28 x 10 -3 in units in which 
c = G = M© = 1. The first column labels the different 
modes. The second column shows the frequencies obtained 
with our code, the third column the results obtained with 
a different nonlinear code based on Cauchy slices [^|. The 
fourth column indicates the frequencies obtained from a lin- 
ear perturbation code ]2^] for the quadrupolar modes. The 
last column shows the relative difference between the present 
code and the results of |23] in percent. 



Mode 


Present Code 


Code (H 


Perturbation 


^| Difference 




(kHz) 


(kHz) 


(kHz) 


(per cent) 


F 


2.690 


2.706 




0.59 


Hi 


4.636 


4.547 




1.96 


Hi 


6.532 


6.320 




3.35 


H 3 


8.418 


8.153 




3.25 


7 


1.388 


1.335 




3.97 




3.504 


3.473 




0.89 




5.510 


5.335 




3.28 


1 P3 


7.400 


7.136 




3.70 


V 


1.871 


1.852 


1.884 


1.03 




4.143 


4.100 


4.110 


1.05 




6.135 


6.019 


6.035 


1.93 


2 P, 


8.087 


7.867 


7.873 


2.80 



using the perturbation 



(87) 



dp = Ap c sin y-^- 

s P = (i+-)p s A 
n p 

where A is the amplitude of the perturbation. 

Additionally, to excite the I = 1,2 non-radial modes 
we perturb the angular velocity component according to 



7TT 

"1 = As[li [ -R2 



u y = Asm [ | < us.' 



(89) 
(90) 



respectively. Fig. |^ shows the initial setup and a contour 
plot of the perturbation given by Eq. (Q) . 

Following [^2|, [2j|, in order to determine the different 
oscillation modes, we analyze the time evolution of dif- 
ferent (fluid and metric) variables at a fixed coordinate 
location. We have checked that the frequencies of the 
oscillation modes are largely independent of the specific 
location. Hence, for the results presented here, we have 
restricted ourselves to extracting those frequencies at co- 
ordinates (x,y) — (¥f-Ax, ^fAy), where N x denotes the 
number of radial zones covering the star. In addition we 
use a radial coordinate 



r = 15- 



1 



(91) 



and we choose a perturbation amplitude of A = 10~ 3 . 
We focus on the time evolution of the radial velocity u x 
for the extraction of radial modes and on the angular 
velocity u y for the extraction of non-radial modes. The 
calculation of mode frequencies follows the work of [p3[ , 
i.e. we determine the zeros of the first derivative of the 
Fourier transform. These zeros correspond to maxima 
in the Fourier transform which are associated with the 
excited modes of oscillation. 



B. Fixed background evolutions 

We start presenting the results for the mode frequen- 
cies of the above stellar model obtained in evolutions in 
which we fix the background geometry of the spacetime 
(i.e. we adopt the so-called Cowling approximation). We 
compare the results to the literature, thus testing the 
hydrodynamics solver of our code. 

Table || shows the results for the fundamental radial 
mode (F) and the first three overtones (Hi — H^) and 
the / and p modes for both, the I = 1 and I = 2 per- 
turbation. Our results are in good agreement with both, 
a different nonlinear code |23fl and an independent linear 
code based upon perturbation theory j22| . There are sev- 
eral reasons which explain the observed differences: First, 
the grid resolution used in our simulations is rather low 
(451 x 21 grid points in x and y, respectively) compared 
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to the finer grids used in p3| (200 x 80), especially in 
the angular direction. For the current setup - in con- 
trast to p3| - only about one half of the radial zones is 
used to cover the star. We note that we use this choice 
here for comparisons with the results of the next section, 
where we also have to resolve gravitational waves in the 
exterior spacetime. Second, we have not implemented 
an atmosphere surrounding the star in a few zones as 
in pH ], which enables the star to radially contract and 
expand. Therefore, the surface of the star may be too 
rigid, which may affect the mode frequencies slightly. 
The numerical implementation of such an atmosphere, 
however, requires a problem-dependent adaptation of its 
parameters. Thi rd, a nd perhaps most importantly, as 
described in Sec. [II A we use a second order reconstruc- 



tion scheme at the cell interfaces, and not the third order 
PPM scheme used in It is worth emphasizing the 

importance of using high-order schemes for the hydrody- 
namics in order to improve the frequency identification 
(see related discussion in [|2|). Nevertheless, for our only 
purpose of assessing the validity of the code we think that 
the overall agreement found is satisfactory. 




2 3 
Bondi time [ms] 



C. Metric-fluid coupled evolutions 

In this section we extract the mode frequencies of the 
above stellar model from fully coupled evolutions of the 
fluid and the geometry. 

Fig. fj] shows the (Bondi) time evolution of the radial 
velocity u x for the I — perturbation (top panel), as well 
as the Fourier transform of this quantity. In the same 
way, Figs. ^ and ^| display the angular velocity evolution 
u v and the Fourier transform for the I = 1 and I = 2 
perturbation, respectively. 

The final evolution time in Figs. corresponds to 5 
ms. The distinctive oscillatory pattern depicted in these 
figures is mainly a superposition of the lowest-order nor- 
mal modes of the fluid. The high-frequency modes are 
usually damped fast by the intrinsic viscosity of the nu- 
merical schemes and at late times the star mostly pul- 
sates in its lowest frequency modes. 

We summarize our results on the mode frequencies in 
Table |J. Note that due to the conservation of linear mo- 
mentum the 1 / mode does not exist (see Fig. ||) (^3|. As 
in the Cowling simulations of the previous section we also 
find now good agreement when comparing to results of 
an independent nonlinear code |39| and to results of lin- 
ear perturbation theory |39], [13 1 . The reasons mentioned 
in the above section for the observed discrepancies are 
still valid here, together with the new source of error in- 
troduced by the metric evolution. 



D. Gravitational waveform 



H1 





H3 



2 4 6 8 10 

frequency [kHz] 

FIG. 7: Pulsations of a K = 100, n = l,p c = 1.28 x 10 -3 
polytrope (c = G — Mq — 1). Top panel: Time evolution 
of the radial velocity u x for the I — perturbation. Bottom 
panel: Fourier transform of the radial velocity. We have la- 
beled the different oscillation modes, F (fundamental) and 
Hi — H3 (first three overtones). 



section concerning the I = 2 perturbation. In Fig. |l0| we 
plot the Bondi news function at the equator N(ys = 0) 
as a function of the observer time. Due to the equatorial 
plane symmetry inherent to the perturbation, which is 
conserved during the evolution, we have vb(jj = 0) = 0, 
which enables us to directly plot N(y = 0), thus avoiding 
suitable interpolations for the wave extraction. 

In order to estimate the amplitude of the signal in 
Fig. [l0| we have also computed the gravitational wave 
emission due to the quadrupole formula. In the light- 
cone approach, the quadrupole news takes the form [031 



To end the validation of our code we study the grav- 
itational wave signal from the simulations of the above 



N = 'Q. 



(92) 
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FIG. 8: Pulsations of a K = 100, n = l,p c = 1.28 x 10" 3 
polytrope (c = G = Mo = 1). Top panel: Time evolution of 
the angular velocity u y for the I = 1 perturbation. Bottom 
panel: Fourier transform of the angular velocity. 



FIG. 9: Pulsations of a K = 100, n = l,p c = 1.28 x 10" 3 
polytrope (c = G = Mo = 1). Top panel: Time evolution of 
the angular velocity u y for the I = 2 perturbation. Bottom 
panel: Fourier transform of the angular velocity. 



The relevant quadrupole moment in axisymmetry reads 

Q = 7:sm 2 e[ dr' [ sm9'd9'r' 4 p(l + e)(^- cos 2 6'- \). 

Jo Jo 2 2 

(93) 

We have included the factor (1 + e) in the formula to ac- 
count for the relativistic internal energy. In the numerical 
calculation of Q we equate the radial and angular coor- 
dinates r' and 9' with the radial and angular coordinates 
r and 9. It is well known [[l5| that the three time deriva- 
tives in Eq. ( |92"|) can cause severe numerical error. In 
order to avoid this problem we fit a cosine to the time 
evolution of the quadrupole moment Q obtained in the 
numerical simulation. The dashed line in Fig. [l^ shows 
the third time derivative of this curve, where the time 
derivatives are taken with repect to Bondi time. 

By evaluating the contributions of the gravitational 
wave signal in Fig. [l^, we see that the dominant contri- 
bution originates from the 2 / mode. The frequency we 



extract from the waveform is 1.57 kHz, in good agreement 
with the result shown in Table ||. There are additional 
contributions to the gravitational wave signal. The exci- 
tation of the l p\ mode is created at the outer boundary 
of the fluid, where we do not allow the star to radially 
contract or to expand. 

Comparing the amplitude of the Bondi news and the 
quadrupole news for the 2 / mode we stress that the am- 
plitude roughly agree. We cannot exclude that the dif- 
ferences are mainly due to numerical errors. Due to the 
different contributions of the time derivative for 7 and 
gauge terms the calculation of N is a difficult task [j30| . 
In addition, we note that the total energy radiated away 
for our setup corresponds only to a tiny fraction of the 
total mass of the spacetime. More precisely, whereas the 
total spacetime mass is IAMq, the total mass radiated 
is only 2.8 x 10~ 9 M Q , which is smaller than the typi- 
cal numerical errors in the determination of the Bondi 
mass. Moreover, most of this radiation is a consequence 
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TABLE II: Mode frequencies obtained in the coupled evo- 
lution for the relativistic polytrope with K = 100, n — 1 
and central density p c = 1.28 x 10 -3 in units in which 
c = G = M a = 1. The first column labels the different modes. 
The second column shows the frequencies obtained with our 
code and the third column shows the results obtained from 
linear perturbation theory Jj9| ^3[ . The last column shows 
the deviations in percent. 



Mode Present Code Perturbation 



Difference 





(kHz) 


(kHz) 


(per cent) 


F 


1.422 


1.442 


1.38 


Hi 


3.993 


3.955 


0.96 


Hz 


6.021 


5.916 


1.77 


Hz 


7.968 


7.776 


2.46 


1 Pi 


1.951 






L p2 


2.844 






1 P, 


5.019 






V 


1.587 


1.579 


0.51 


V 


3.748 


3.710 


1.02 


2 P2 


5.811 


5.689 


2.14 


2 P3 


7.848 


7.580 


3.54 



Bondi news 
averaged Bondi news 
quadrupole formula 




does not take into account the curvature of spacetime. 

At late times (after about 3 ms), the Bondi news does 
not strictly oscillate around zero. This numerical ef- 
fect can be weakened using more sophisticated numerical 
methods for the fluid update (e.g. third-order cell recon- 
struction procedures). We note that similar drifts are re- 
ported in |2^, |39) in time evolutions of equilibrium models 
using a perfect fluid EoS. For the results shown in Fig. [l(], 
we have used a polytropic EoS during the fluid evolution. 
This is legitimate as the effect of heating is negligible for 
our stellar object close to equilibrium. Finally, we note 
that the maximum deviation from global energy conser- 
vation, Eq. (p2|), in this simulation is 2.5x 10~ 5 M Q , which 
is very adequate for the long integration time. 



VI. CONCLUSION 

We have presented numerical algorithms to solve the 
coupled Einstein-perfect fluid system in axisymmctry. 
Our approach is based upon the characteristic formula- 
tion of general relativity in which spacetime is foliated 
with a family of outgoing light cones emanating from a 
regular center. Due to a suitable compactification of the 
spacetime future null infinity is part of our finite numer- 
ical grid where we unambiguously extract gravitational 
waves. 

Applying the nonlinear, fully relativistic code to stud- 
ies of neutron stars modeled as polytropes, it has passed 
several tests, aimed at testing both the fluid evolution as 
well as the metric solver in the nonlinear regime. The 
code can accurately maintain long-term stability of rel- 
ativistic stars and we have applied it to the study of 
stellar pulsations. We have extracted the frequencies 
of different non-radial fluid modes and the gravitational 
wave signal arising in time evolutions of perturbed stellar 
configurations. Applications of the present code in the 
computation of the gravitational waveforms emitted in 
axisymmetric core collapse situations will be presented 
elsewhere. 



Bondi time [ms] 



FIG. 10: News function as a function of observer time at in- 
finity at the equator. The (unphysical) initial signal is due 
to the gravitational wave content in the initial data. The 
main oscillation frequency corresponds to the 2 / mode with 
other frequencies overlayed. The dotted curve shows a suit- 
able average of the numerical Bondi news smearing out higher 
frequencies. The dashed line corresponds to the estimated 
amplitude of the 2 / mode oscillation which is extracted using 
the quadrupole formula. See text for a detailed discussion. 

of the initial gravitational wave content and is radiated 
away in between ub — ms and ub = 0.1 ms. Addi- 
tional differences in the curves for the Bondi news and 
the quadrupole news in Fig. [l^ might be created by the 
rough estimate of the quadrupole Q in Eq. ( p3|) which 
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